%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Boucle et post traitement du poblème de fissure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% ATTENTION !!!
% Il faut rentrer les valeurs de H1, LF et YCF avant de lancer les calculs !!!!



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

setenv('CASTEM','/home/andres/Cast3m/' ) ;


H1=10  ;
H2=10  ;
YCF=0  ;
LF=10  ;
L=50   ;

% H1=10 /1000;
% H2=10 /1000;
% YCF=0 /1000;
% LF=10 /1000;
% L=50 /1000;

materiau.E = 300e9 ;
materiau.nu = 0.3 ;
materiau.lambda = materiau.E*materiau.nu/((1+materiau.nu)*(1-2*materiau.nu)) ;
materiau.mu = materiau.E/(2*(1+materiau.nu)) ;

geometrie.longueur=L;
geometrie.hauteur=H1+H2;


angleInf = 2;
angleMax = 2;
pas = 5;

angleCalcul = angleInf:pas:angleMax;

theta_calcule=zeros(size(angleCalcul,1),1);

DELIMITER = ';';
HEADERLINES = 0;

for i = 1: length(angleCalcul)

    fid = fopen('Fissure.dgibi','r+');
    fprintf(fid, 'ANGLEF = %2.2f ; \n', angleCalcul(i));
    
    system('castem11 Fissure.dgibi');
    
    
 %---------------------------------------------------------------------   
    calcul.x_haut = importdata('x_haut', DELIMITER, HEADERLINES) ; %/1000;
    
    calcul.x_gauche1 = importdata('x_gauche1', DELIMITER, HEADERLINES) ; %/1000;
    
    calcul.x_gauche2 = importdata('x_gauche2', DELIMITER, HEADERLINES) ; %/1000;
    
    calcul.x_droite1 = importdata('x_droite1', DELIMITER, HEADERLINES) ; %/1000;
    
    calcul.x_droite2 = importdata('x_droite2', DELIMITER, HEADERLINES) ; %/1000;
    
    %.................................................................
    
    calcul.y_haut = importdata('y_haut', DELIMITER, HEADERLINES) ; %/1000;
    
    calcul.y_gauche1 = importdata('y_gauche1', DELIMITER, HEADERLINES) ; %/1000;
    
    calcul.y_gauche2 = importdata('y_gauche2', DELIMITER, HEADERLINES) ; %/1000;
    
    calcul.y_droite1 = importdata('y_droite1', DELIMITER, HEADERLINES) ; %/1000;
    
    calcul.y_droite2 = importdata('y_droite2', DELIMITER, HEADERLINES) ; %/1000;
    
 %--------------------------------------------------------------------- 
    calcul.sigmaxx_haut = importdata('sigmaxx_haut', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmaxx_gauche1 = importdata('sigmaxx_gauche1', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmaxx_gauche2 = importdata('sigmaxx_gauche2', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmaxx_droite1 = importdata('sigmaxx_droite1', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmaxx_droite2 = importdata('sigmaxx_droite2', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmaxx_bas = importdata('sigmaxx_bas', DELIMITER, HEADERLINES) ; %*1e6;
    
    %.................................................................
    
    calcul.sigmayy_haut = importdata('sigmayy_haut', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmayy_gauche1 = importdata('sigmayy_gauche1', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmayy_gauche2 = importdata('sigmayy_gauche2', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmayy_droite1 = importdata('sigmayy_droite1', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmayy_droite2 = importdata('sigmayy_droite2', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmayy_bas = importdata('sigmayy_bas', DELIMITER, HEADERLINES) ; %*1e6;
    
    %.................................................................
    
    calcul.sigmaxy_haut = importdata('sigmaxy_haut', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmaxy_gauche1 = importdata('sigmaxy_gauche1', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmaxy_gauche2 = importdata('sigmaxy_gauche2', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmaxy_droite1 = importdata('sigmaxy_droite1', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmaxy_droite2 = importdata('sigmaxy_droite2', DELIMITER, HEADERLINES) ; %*1e6;
    
    calcul.sigmaxy_bas = importdata('sigmaxy_bas', DELIMITER, HEADERLINES) ; %*1e6;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    resultats(i)=calcul;
    
    
    % Adadptation des données
    
%     mihauteur_gauche=H1+YCF-(LF*sin(i *pi/180));
%     mihauteur_droite=H1+YCF+(LF*sin(i *pi/180));


    
    % haut : OK
    % bas : OK
    % droite1 : OK
    % gauche1 : OK
    
%     % droite2
%     
%     calcul.x_droite2(:,1) = calcul.x_droite2(:,1) + mihauteur_droite ;
%     calcul.y_droite2(:,1) = calcul.y_droite2(:,1) + mihauteur_droite ;
%     calcul.sigmaxx_droite2(:,1) = calcul.sigmaxx_droite2(:,1) +  mihauteur_droite ;
%     calcul.sigmaxy_droite2(:,1) = calcul.sigmaxy_droite2(:,1) +  mihauteur_droite ;
%     calcul.sigmayy_droite2(:,1) = calcul.sigmayy_droite2(:,1) +  mihauteur_droite ;
%     
%     % gauche2
%     
%     calcul.x_gauche2(:,1) = calcul.x_gauche2(:,1) + mihauteur_gauche ;
%     calcul.y_gauche2(:,1) = calcul.y_gauche2(:,1) + mihauteur_gauche ;
%     calcul.sigmaxx_gauche2(:,1) = calcul.sigmaxx_gauche2(:,1) +  mihauteur_gauche ;
%     calcul.sigmaxy_gauche2(:,1) = calcul.sigmaxy_gauche2(:,1) +  mihauteur_gauche ;
%     calcul.sigmayy_gauche2(:,1) = calcul.sigmayy_gauche2(:,1) +  mihauteur_gauche ;
    
    % Extraction des données
    
    x_haut = calcul.x_haut;
    x_gauche1 = calcul.x_gauche1;
    x_gauche2 = calcul.x_gauche2;
    x_droite1 = calcul.x_droite1;
    x_droite2 = calcul.x_droite2;
    
    y_haut = calcul.y_haut;
    y_gauche1 = calcul.y_gauche1;
    y_gauche2 = calcul.y_gauche2;
    y_droite1 = calcul.y_droite1;
    y_droite2 = calcul.y_droite2;
    
    sigmaxx_haut = calcul.sigmaxx_haut;
    sigmaxx_bas = calcul.sigmaxx_bas;
    sigmaxx_gauche1 = calcul.sigmaxx_gauche1;
    sigmaxx_gauche2 = calcul.sigmaxx_gauche2;
    sigmaxx_bas(:,2) = -calcul.sigmaxx_bas(:,2);  %%
    sigmaxx_gauche1(:,2) = -calcul.sigmaxx_gauche1(:,2);  %%
    sigmaxx_gauche2(:,2) = -calcul.sigmaxx_gauche2(:,2);  %%
    sigmaxx_droite1 = calcul.sigmaxx_droite1;
    sigmaxx_droite2 = calcul.sigmaxx_droite2;
    
    sigmaxy_haut = calcul.sigmaxy_haut;
    sigmaxy_bas = calcul.sigmaxx_bas;
    sigmaxy_gauche1 = calcul.sigmaxy_gauche1;
    sigmaxy_gauche2 = calcul.sigmaxy_gauche2;
    sigmaxy_bas(:,2) = -calcul.sigmaxx_bas(:,2);  %%
    sigmaxy_gauche1(:,2) = -calcul.sigmaxy_gauche1(:,2);  %%
    sigmaxy_gauche2(:,2) = -calcul.sigmaxy_gauche2(:,2);  %%
    sigmaxy_droite1 = calcul.sigmaxy_droite1;
    sigmaxy_droite2 = calcul.sigmaxy_droite2;
    
    sigmayy_haut = calcul.sigmayy_haut;
    sigmayy_bas = calcul.sigmayy_bas;
    sigmayy_gauche1 = calcul.sigmayy_gauche1;
    sigmayy_gauche2 = calcul.sigmayy_gauche2;
    sigmayy_bas(:,2) = -calcul.sigmayy_bas(:,2);  %%
    sigmayy_gauche1(:,2) = -calcul.sigmayy_gauche1(:,2);  %%
    sigmayy_gauche2(:,2) = -calcul.sigmayy_gauche2(:,2);  %%
    sigmayy_droite1 = calcul.sigmayy_droite1;
    sigmayy_droite2 = calcul.sigmayy_droite2;
    
    
    
    mihauteur_gauche=x_gauche1(end,1);
    mihauteur_droite=x_droite1(end,1);
    
    
    % droite2
    
    x_droite2(:,1) = x_droite2(:,1) + mihauteur_droite ;
    y_droite2(:,1) =  y_droite2(:,1) + mihauteur_droite ;
    sigmaxx_droite2(:,1) = sigmaxx_droite2(:,1) +  mihauteur_droite ;
    sigmaxy_droite2(:,1) = sigmaxy_droite2(:,1) +  mihauteur_droite ;
    sigmayy_droite2(:,1) = sigmayy_droite2(:,1) +  mihauteur_droite ;
    
    % gauche2
    
    x_gauche2(:,1) = x_gauche2(:,1) + mihauteur_gauche ;
    y_gauche2(:,1) = y_gauche2(:,1) + mihauteur_gauche ;
    sigmaxx_gauche2(:,1) = sigmaxx_gauche2(:,1) +  mihauteur_gauche ;
    sigmaxy_gauche2(:,1) = sigmaxy_gauche2(:,1) +  mihauteur_gauche ;
    sigmayy_gauche2(:,1) = sigmayy_gauche2(:,1) +  mihauteur_gauche ;
    
    % Calcul de l'écart à la réciprocité
    
%         donnees.geom=geometrie;
%     donnees.mat=materiau;
%     donnees.calcul=calcul;
%     
%     Ex = ER(donnees,'x') ;
%     Ey = ER(donnees,'y') ;
%     Exy = ER(donnees,'xy') ;
%     
% Ex=ER( materiau.lambda,materiau.mu,geometrie.hauteur,geometrie.longueur,'x',x_haut,x_gauche1,x_gauche2,...
%     x_droite1,x_droite2,y_haut,y_gauche1,y_gauche2,y_droite1,y_droite2,...
%     sigmaxx_haut,sigmaxx_bas,sigmaxx_gauche1,sigmaxx_gauche2,sigmaxx_droite1,...
%     sigmaxx_droite2,sigmaxy_haut,sigmaxy_bas,sigmaxy_gauche1,sigmaxy_gauche2,...
%     sigmaxy_droite1,sigmaxy_droite2,sigmayy_haut,sigmayy_bas,sigmayy_gauche1,...
%     sigmayy_gauche2,sigmayy_droite1,sigmayy_droite2);
% 
% Ey=ER( materiau.lambda,materiau.mu,geometrie.hauteur,geometrie.longueur,'y',x_haut,x_gauche1,x_gauche2,...
%     x_droite1,x_droite2,y_haut,y_gauche1,y_gauche2,y_droite1,y_droite2,...
%     sigmaxx_haut,sigmaxx_bas,sigmaxx_gauche1,sigmaxx_gauche2,sigmaxx_droite1,...
%     sigmaxx_droite2,sigmaxy_haut,sigmaxy_bas,sigmaxy_gauche1,sigmaxy_gauche2,...
%     sigmaxy_droite1,sigmaxy_droite2,sigmayy_haut,sigmayy_bas,sigmayy_gauche1,...
%     sigmayy_gauche2,sigmayy_droite1,sigmayy_droite2);
% 
% Exy=ER( materiau.lambda,materiau.mu,geometrie.hauteur,geometrie.longueur,'xy',x_haut,x_gauche1,x_gauche2,...
%     x_droite1,x_droite2,y_haut,y_gauche1,y_gauche2,y_droite1,y_droite2,...
%     sigmaxx_haut,sigmaxx_bas,sigmaxx_gauche1,sigmaxx_gauche2,sigmaxx_droite1,...
%     sigmaxx_droite2,sigmaxy_haut,sigmaxy_bas,sigmaxy_gauche1,sigmaxy_gauche2,...
%     sigmaxy_droite1,sigmaxy_droite2,sigmayy_haut,sigmayy_bas,sigmayy_gauche1,...
%     sigmayy_gauche2,sigmayy_droite1,sigmayy_droite2);
% 
%     
% 
%     % Exploitation des résultats
%     
%     Un = -(Ex+Ey)/(2*(materiau.lambda + materiau.mu)) ;  % il manque un moins
%     Unn = (Ey-Ex)/(2*materiau.mu) ;
%     Utt = -Exy/materiau.mu ;
%     U = sqrt( Unn^2 + Utt^2 ) ;
%     Ut = sqrt( U^2 - Un^2 ) ;
%     
%     [ cs_alpha ] = [ Un , -Ut ; Ut , Un ] \ [ Unn ; Utt ] ;
%     
%     c_alpha=cs_alpha(1);
%     s_alpha=cs_alpha(2);
%     alpha_calcule = (180/pi)*acos(c_alpha) ;
% %     alpha_calcule = (180/pi)*asin(s_alpha) ;
%     
%     theta_calcule(i) = -alpha_calcule/2 + 90  % angle de la fissure, pas de la normale
%     

E11=ER( materiau.E,materiau.nu,geometrie.hauteur,geometrie.longueur,'Andrieux_11',x_haut,x_gauche1,x_gauche2,...
    x_droite1,x_droite2,y_haut,y_gauche1,y_gauche2,y_droite1,y_droite2,...
    sigmaxx_haut,sigmaxx_bas,sigmaxx_gauche1,sigmaxx_gauche2,sigmaxx_droite1,...
    sigmaxx_droite2,sigmaxy_haut,sigmaxy_bas,sigmaxy_gauche1,sigmaxy_gauche2,...
    sigmaxy_droite1,sigmaxy_droite2,sigmayy_haut,sigmayy_bas,sigmayy_gauche1,...
    sigmayy_gauche2,sigmayy_droite1,sigmayy_droite2);

E22=ER( materiau.E,materiau.nu,geometrie.hauteur,geometrie.longueur,'Andrieux_22',x_haut,x_gauche1,x_gauche2,...
    x_droite1,x_droite2,y_haut,y_gauche1,y_gauche2,y_droite1,y_droite2,...
    sigmaxx_haut,sigmaxx_bas,sigmaxx_gauche1,sigmaxx_gauche2,sigmaxx_droite1,...
    sigmaxx_droite2,sigmaxy_haut,sigmaxy_bas,sigmaxy_gauche1,sigmaxy_gauche2,...
    sigmaxy_droite1,sigmaxy_droite2,sigmayy_haut,sigmayy_bas,sigmayy_gauche1,...
    sigmayy_gauche2,sigmayy_droite1,sigmayy_droite2);

E12=ER( materiau.E,materiau.nu,geometrie.hauteur,geometrie.longueur,'Andrieux_12',x_haut,x_gauche1,x_gauche2,...
    x_droite1,x_droite2,y_haut,y_gauche1,y_gauche2,y_droite1,y_droite2,...
    sigmaxx_haut,sigmaxx_bas,sigmaxx_gauche1,sigmaxx_gauche2,sigmaxx_droite1,...
    sigmaxx_droite2,sigmaxy_haut,sigmaxy_bas,sigmaxy_gauche1,sigmaxy_gauche2,...
    sigmaxy_droite1,sigmaxy_droite2,sigmayy_haut,sigmayy_bas,sigmayy_gauche1,...
    sigmayy_gauche2,sigmayy_droite1,sigmayy_droite2);



    % Exploitation des résultats

%     Ka = -(E11+E22)
%     Kb = -(E22-E11)
%     Kc = -2*E12
%     
%     theta_calcule(i)=intersec(Ka,Kb,Kc);
    
U=sqrt((E11-E22)^2+E12^2);
theta=acos(sqrt(U^2-(E11+E22)^2)/U)*180/pi



%     pause
    
    
    
    
    
end

% theta_calcule
